function output = my_getLongTermDetrended(Data,Time,opt)
    % output = getLongTermDetrended(Data)
    % 函数说明：使用最小二乘线性拟合，去掉数据中的长期趋势
    % 函数输入：
    %   Data：二维或三维矩阵，二维矩阵需要时间序列的列向量构成，三维矩阵需要第三维是
    %       时间维度。
	%   Time: option, needed if remove montly trend
	%   option:
	%   	use as : 'key',value
	%   	'MonthlyLongTerm': true or false, trigger for remove long term trend in 
	%   		each month
	%		'MonthlyPeriod': true or false, trigger for remove seasonal trend
	%		'nanflag': 'includnan' or 'omitnan', including nan or not
    % log:
	% 	2022/03/21 update support nanflag
	% 	2022/03/25 ver-0.3 support monthly detrend
	arguments
		Data
		Time = NaT
		opt.nanflag = 'includenan'
		opt.MonthlyLongTerm = false
		opt.MonthlyPeriod = false
	end

	% load data
	nanflag = opt.nanflag;

    switch ndims(Data)
    case 2
        output = detrend(Data,nanflag);
		if opt.MonthlyLongTerm || opt.MonthlyPeriod
			output = core_monthly(Data,Time,...
				nanflag,...
				opt.MonthlyLongTerm,...
				opt.MonthlyPeriod);
		end
    case 3
        [Data2D, SIZE] = my_shiftdim(Data, 1);
        output2D = detrend(Data2D,1,nanflag);
		if opt.MonthlyLongTerm || opt.MonthlyPeriod
			output2D = core_monthly(output2D,Time,...
				nanflag,...
				opt.MonthlyLongTerm,...
				opt.MonthlyPeriod);
		end
        output = my_shiftdim(output2D,2,SIZE);
    end
end

function out =  core_monthly(Data,Time,nanflag,MonthlyLongTerm,MonthlyPeriod)
	% detrend:
	% 	MonthlyLongTerm: each month has long term trend, trend is different
	% 		between months
	% 	MonthlyPeriod: seasonal trend, removed by minus average mean of each month
	% Input:
	% 	Data: 2D matrix, 2nd Dimension is time dimension
	% 	nanflag: whether include nan while using function detrend
	% 	MonthlyLongTerm: logical, trigger for detrend
	% 	MonthlyPeriod: logical, trigger for detrend
	% output:
	% 	out: data matrix(time-site matrix) removed trend
	
	out = NaN(size(Data));
	for ii = 1:size(Data,2)
		ts = Data(:,ii);
		% reshape month, to calculate per month
		[ts_reshaped,Indice] = reshapeByMonth(ts,Time);

		% remove long term trend
		if MonthlyLongTerm
			% 1 is n degree trend
			ts_reshaped = detrend(ts_reshaped,1,nanflag);
		end

		% remove seasonal trend
		if MonthlyPeriod
			ts_reshaped = removeMonthlyTrend(ts_reshaped);
		end

		% reshape month, return to mornal format
		out(:,ii) = ts_reshaped(Indice);
	end
end

function output = removeMonthlyTrend(TsReshaped)
	% remove long term trend each month
	% output is still reshaped month, col is month,row is year
	TrendBase = mean(TsReshaped,1,'omitnan');
	Trend = repmat(TrendBase,size(TsReshaped,1),1);
	output = TsReshaped - Trend;
end

function [output,Indice] = reshapeByMonth(input,Time)
	% 将单一变量的时间序列，按照月份排序，组成列为同一月份,行为同一年的矩阵
	% 	2001,1 2001,2 2001,3 ... 2001,12
	% 	2002,1 2002,2 2002,3 ... 2002,12
	% 	...
	% input是时间序列
	% Time是datetime格式，跟input是相同的大小
	% 不足的用nan补全，所以后续计算的时候需要omitnan
	% 要求序列本身没有nan值，不然无法区分
	
	% 检查开头的月份
	checkSize(input);
	input = input(:);
	Head = ones(Time(1).Month - 1,1) + nan;
	Tail = ones(12 - Time(end).Month,1) + nan;
	output = [Head;input(:);Tail];
	output = reshape(output,12,[]);
	output = output';

	Indice = 1:numel(output);
	Indice(1:numel(Head)) = nan;
	Indice((end-numel(Tail)+1):end) = nan;
	Indice = reshape(Indice,size(output));

	Indice = Indice';
	Indice = Indice(:);
	Indice = Indice(~isnan(Indice));


	function checkSize(input)
		% check input is vector of time series
		if isvector(input)
			true;
		else
			error('input of reshapeByMonth must be vector');
		end
	end
end
